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In conjunction with NASA and the Department of Defense, the Johns Hopkins Univer- 
sity Applied Physics Laboratory (JHU/APL) has been investigating analytical techniques 
to address many of the fundamental issues associated with solar exploration spacecraft and 
high-speed atmospheric vehicle systems. These issues include: thermo-structural response 
including the effects of thermal management via the use of surface optical properties for 
high-temperature composite structures; aerodynamics with the effects of non-equilibrium 
chemistry and gas radiation; and aero-thermodynamics with the effects of material ablation 
for a wide range of thermal protection system (TPS) materials. The need exists to integrate 
these discrete tools into a common framework that enables the investigation of interdisci- 
plinary interactions (including analysis tool, applied load, and environment uncertainties) 
to provide high fidelity solutions. 

In addition to developing robust tools for the coupling of aerodynamically induced ther- 
mal and mechanical loads, JHU/APL has been studying the optimal design of high-speed 
vehicles as a function of their trajectory. Under traditional design methodology the opti- 
mization of system level mission parameters such as range and time of flight is performed 
independently of the optimization for thermal and mechanical constraints such as stress 
and temperature. A truly optimal trajectory should optimize over the entire range of 
mission and thermo-mechanical constraints. 

Under this research, a framework for the robust analysis of high-speed spacecraft and 
atmospheric vehicle systems has been developed. It has been built around a generic, loosely 
coupled framework such that a variety of readily available analysis tools can be used. The 
methodology immediately addresses many of the current analysis inadequacies and allows 
for future extension in order to handle more complex problems. 

I. Introduction 

In order to perform a full transient aero-thermo-mechanical trajectory analysis three different physics- 
based models must be coupled: a six degree-of-freedom (6D0F) kinematics model, a computational fluid 
dynamics (CFD) model, and a computational structural mechanics (CSM) model. A range of coupling 
technologies are available and must be evaluated on a per- application basis. Direct coupling involves solving 
a complete, fully coupled set of equations that simultaneously solve for all of the necessary fluid and structural 
field variables. This type of coupling is computationally expensive since the most limiting time and length 
scales will drive the simulation and the equations must be solved over a very large domain which encompasses 
both the vehicle and the surrounding fluid. There are currently no commercial applications of this type so 
solutions using this method will require substantial code development and be limited, at least in the near- 
term, by a lack of robustness. An alternative to co-solution is a server-client architecture where the 6D0F, 
CFD, and CSM codes are linked through a central server which manages the translation of data across the 
simulation boundaries. Mesh-based Parallel Code Coupling Interface (MpCCI) developed by Fraunhofer- 
Institute SCAI is a widely supported server-client coupling tool which advertises an open architecture; 
however, the server and the application programming interface (API) are complex and necessitate a high 
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degree of mastery of not only each coupled code but also of the server tool itself. Recently, several of the 
commercial codes that actively supported the MpCCI interface have begun developing their own tightly 
coupled CFD-CSM capabilities. A third option is loose coupling which allows interpolation of data from a 
CFD solution to a CSM application based upon a trajectory computed from a 6DOF simulation. In this 
scheme the individual codes are coupled but in a uni-directional fashion. Thus, fluid dynamic pressures and 
fluxes depend upon the CSM surface temperatures and the 6DOF generated flight conditions, but geometry 
and mass property changes from the thermo-mechanical solution do not affect the CFD or 6DOF solutions. 
For a given application, the degree of coupling necessary to obtain solutions of a desired fidelity needs to be 
evaluated; however, for the applications of interest in this research loose coupling is believed to be reasonably 
accurate and substantially more efficient. Thus, a loose coupling technique has been employed. A functional 
diagram of the loose coupling procedure is shown in Figure 1 



Figure 1. Functional diagram of the loose coupling process. 


where the tasks required to perform the optimized trajectory analysis are labeled sequentially and described 
below. Currently, this methodology is restricted to using the user subroutines available in ABAQUS; however, 
it could easily be extended to other CSM codes as long as they have similar functionality. 

1. CAD: Based upon the current or desired design a computer aided design (CAD) model must be created. 
The same CAD model should be used for both the CFD and CSM calculations. As will be explained in 
following section, the CFD and CSM meshes do not need to be continuous along the shared boundary, 
but in order to insure that the forces are transferred properly, it is extremely important that the solid 
model boundaries coincide. 

2. CFD: Because this method is loosely coupled it is not necessary to recompute the CFD results during 
the simulation. Thus, the CFD data can be pre-computed and stored in a database that is queried as 
a function of the current vehicle flight conditions. For this work, the database may be a function of the 
altitude, total angle- of- attack, aerodynamic roll angle, and Mach number. The total angle-of-attack 
and aerodynamic roll angles are described in Figure 2. where u , v, and w are the vehicle body axes. 



Figure 2. Aerodynamic Angle Definitions. 


3. Optimizer: The 6DOF and the CSM codes are coupled within the optimizer step. The 6DOF reads 
force and moment coefficients from the pre-computed CFD analysis and then based upon a set of initial 
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conditions propagates the vehicle in time and space until a predetermined time or position is reached. 
The CSM code reads the time history data from the 6DOF and using the CFD database determines the 
appropriate mechanical (pressure) and thermal (flux) boundary conditions as a function of time along 
the trajectory. Prior to launching the simulation optimization variables and corresponding thresholds 
are established. After the CSM analysis is complete for a given trajectory the thresholds are checked 
by the optimizer and the initial conditions are iterated until a converged solution is achieved. 

4. Optimal Trajectories: Ideally, after a small number of iterations a trajectory optimized on a set of 
predetermined conditions will be found. 

The results in this paper were obtained using CFD++ as the CFD code and ABAQUS as the CSM code. 

The 6DOF and the Optimizer were developed in-house. 


II. Mesh Interpolation 


The CFD and CSM solvers need not and, in most cases, should not have matching boundary grids along 
the fluid-structure interface. It is, therefore, necessary to develop an interpolation scheme to translate the 
CFD data to the CSM grid. The concept of a kernel approximation for randomly distributed particles that 
has been actively developed in support of the meshless Galerkin and particle-based computational methods 
is used to perform the interpolation. According to this concept, a corrected kernel approximation / (x) to a 
variable / ( x ) at a point x is given by (1) 


f, x_ f 

f[Xi> f w (xi) dn 


(i) 


If a particle approximation of Equation (1) is formed then the new particle-based corrected kernel ap- 
proximation of / (x) is, 


/ ( x i) = 


E * =1 f(x?)W a ( Xj)AW 
£a=l W a (x i) AW 


( 2 ) 


where subscript a identifies the particle at xf on the surface for which / (xf) is known and Xi is the point 
where the interpolated value is desired. In this analysis, xf is a point on the CFD boundary and Xi is an 
integration point on the CSM boundary. W a (x^) is the value of the kernel or weight from particle a at 
point a^and A is a volume associated with particle a. Since the CFD boundary particles are not associated 
with any structural volume, AQ a is assumed to be constant for all a and, thus, (2) becomes the standard 
weighted average given below. 


Eq=l f(*?)W a (Xj, 


A variety of forms of the kernel function are possible and a piecewise 


cubic function was chosen 2 


(3) 


{ (2 — ip) 3 — 4 (1 — ip) 3 0 < ip < 1 
(2-^) 3 1 < ip <2 (4) 

0 ip > 2 

where ip = — ~ Xx ^ . h is the kernel smoothing length and determines the extent of the influence of a 
given particle a. The method of computing h is largely arbitrary; however, it should not be so large that 
the interpolation loses its local character nor should it be so small that there are an insufficient number of 
particles to provide an accurate interpolation. For this analysis, h was chosen such that the influence domain 
of a given particle was a sphere with radius equal to 1.5 times the distance to the nearest neighbor particle. 1 
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III. Time Interpolation 


Prior to starting the analysis a CFD database of pressures (p), temperatures (#), and fluxes ( q ) along the 
boundary must be computed at discrete values of altitude ( h ), total angle-of-attack (<a), aerodynamic roll 
angle (</>), and Mach number (M) for both an adiabatic and a cold- wall condition along the fluid-structural 
boundary. The cold-wall temperatures specified in the cold-wall CFD analysis should be reasonable and 
in the same range expected from the trajectory simulation. For the cases studied in this report sea-level 
ambient temperature provided good results. At the beginning and end of each computation step a three 
variable interpolation is performed as a function of a, 0, and M using an eight node three-dimensional 
interpolation. A second interpolation is then performed to find the data at the correct altitude, h. During 
the CSM analysis step the current value of p, 0 , and q will be linearly interpolated using the values from the 
beginning and the end of the step. More detail on this process is provided in the following derivation. 

At any point on the trajectory it is assumed that the proper conditions may be approximated by the 
interpolation, 


= N}(a,4,M)Zi (5) 

e{*AM) = ( 6 ) 

where the superscripts indicate the nearest CFD database altitude above and below the current actual 
altitude, £ is a generic variable that represents any of the database values (i.e. h , #, 0), and repeated indices 
imply summation. The shape functions, N i: are determined using the standard tri-linear approximation 
given in (7). 


£ {oti 0, — clq T cl\(x T $20 T T < 24 otcf) T CL^ctJVf T < 2 g£)Af T < 27 (y,(f)Ad (7) 

A linear interpolation between the two database altitudes as given in Equation (8) produces the the 
proper interpolated value for £ at the given flight conditions. 

<-{h,e,e)=se+{i-s)e ( 8 ) 

S is a non-dimensional measure of the current altitude as given by, 

S = (h — h 1 ) / ( h 2 — h 1 ) . (9) 

A final linear interpolation between the two time steps, £ n and t n + 1 , then determines the correct values of p, 
#, and q at any time over the time step. A graphical representation of the CSM analysis process is shown in 
Figure 3. 



Time Integration Loop 


Figure 3. CSM Time Integration Procedure Functional Diagram. 


IV. Aero-thermal Heating 

Convective heat transfer through the boundary layer is calculated using the pre-computed CFD adiabatic 
and isothermal cold- wall data. The convective heating relation is estimated using the wall temperature from 
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the adiabatic analysis, O aw , as the driving force. 3 This, when combined with the heat flux from the isothermal 
cold- wall analysis, provides the convection coefficient, /i, as in (10). 


» = 


Qw 


0 W - 0 


aw 


( 10 ) 


Subscript w indicates values taken at the boundary (wall) and subscript a indicates results derived from the 
adiabatic analysis. The convective heat flux is then found through (11). 


Q convective ( @aw @w ) (ii) 

The radiative heat flux may also be applied using the standard radiative heat equation given in (12). 

Qradiative = <?£ “ &t) (12) 

In (10) and (12), /i, q w , and 0 aw are derived from the CFD results, 0 W is found from the CSM solution, 
and <9oo is a function of the trajectory as determined by the 6DOF. It is also important to note that both the 
ambient wall temperature and the heating coefficient are interpolated to the CSM integration points using 
the previously described time accurate interpolation. 


V. Test Case 

The analysis methodology was applied to a simple axisymmetric ballistic re-entry vehicle. A schematic 
of the vehicle is given in Figure 4. 

^ L ^ 

D 

Figure 4. Axisymmetric example case geometry (not to scale). 



The vehicle is solid tungsten with the following material properties: 

p = 19, 300 kg/m 3 
c p = 132.0J/kgK 
k = 174.0 W/mK 
£ = 0.7 

where p is the mass density, c p is the specific heat, k is the thermal conductivity, and £ is the emissivity. 
The total angle-of- attack and the aerodynamic roll angles were zero throughout the trajectory so that the 
CFD database was solely a function of altitude and Mach number. 

V.A. Aerothermal Heating 

A complete ballistic trajectory analysis was performed and compared to results using heating correlations. 4 
For this simple geometry the correlations should be accurate and are, therefore, considered a reasonable 
source for validation. The trajectory altitude and Mach number as a function of time are given in Figure 5. 

The temperatures and stagnation point heat flux are given in Figures 6 and 7, respectively, where the 
data identified as “Lear, Green” was generated using the current methodology and “Detra, Kemp, Riddell” 
was generated using the correlations. The discrete jumps in the heat flux are a result of assuming sea- level 
ambient conditions at all altitudes in the CFD database. Even with this approximation, the results show an 
excellent agreement with the heating correlations. 
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Figure 5. Example Case Trajectory Altitude and Mach Number. 



Figure 6. Stagnation Point Heat Flux. 



Figure T. Stagnation Point Temperature. 
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V.B. Optimization 


A simplified, single variable optimization study was conducted. For this analysis a zero angle-of- attack 
trajectory neglecting gravitational accelerations and using sea-level flight conditions was performed. A 
desired stagnation point temperature was specified and initial velocity was varied using a simple Newton- 
Raphson scheme. The final optimal initial velocity was found to be less than 85% of the initial velocity. The 
analysis converged in 5 iterations and a plot of the optimization process is given in Figure 8. 



Figure 8. Optimization Progress Chart for the Example Case. 


VI. Conclusion 

This research demonstrates that a robust loosely coupled aero-thermomechanical analysis methodology 
can be used with existing commercial codes. This allows the most flexibility in the design process because it 
exploits the robustness of existing codes and the availability of existing code-specific expertise. Additionally, 
this research has shown that this tool can be wrapped within an optimization layer in order to develop vehicle 
trajectories that are optimized on a complete set of variables whether they be performance or structurally 
based. Because of the modular nature of the methodology it is easily extendable to include other physical 
effects and to solve more complex problems. 
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